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Abstract — The problem of clock offset estimation in a two 
way timing message exchange regime is considered when the 
likelihood function of the observation time stamps is Gaussian, 
exponential or log-normally distributed. A parametrized solution 
to the maximum likelihood (ML) estimation of clock offset, based 
on convex optimization, is presented, which differs from the 
earlier approaches where the likelihood function is maximized 
graphically. In order to capture the imperfections in node 
oscillators, which may render a time-varying nature to the clock 
offset, a novel Bayesian approach to the clock offset estimation is 
proposed by using a factor graph representation of the posterior 
density. Message passing using the max-product algorithm yields 
a closed form expression for the Bayesian inference problem. 
Several lower bounds on the variance of an estimator are 
derived for arbitrary exponential family distributed likelihood 
functions which, while serving as stepping stones to benchmark 
the performance of the proposed clock offset estimators, can be 
useful in their own right in classical as well Bayesian param- 
eter estimation theory. To corroborate the theoretical findings, 
extensive simulation results are discussed for classical as well 
as Bayesian estimators in various scenarios. It is observed that 
the performance of the proposed estimators is fairly close to the 
fundamental limits established by the lower bounds. 

Index Terms — Clock synchronization, factor graphs, message 
passing, estimation bounds, wireless sensor networks. 

I. Introduction 

WIRELESS sensor networks (WSNs) typically consist 
of a large number of geographically distributed sensor 
nodes, deployed to observe some phenomenon of interest. The 
nodes constituting such a network are low cost sensors that 
have limited abilities of data processing and communication. 
WSNs envisage tremendous applications in such diverse areas 
as industrial process control, battlefield surveillance, health 
monitoring, target localization and tracking, etc., With 
the recent advances in digital circuit technology, WSNs are 
expected to play a pivotal role in future wireless communica- 
tions. 

Clock synchronization in sensor networks is a critical 
component in data fusion and duty cycling operations, and 
has gained widespread interest in the past few years. Most 
of the current methods consider sensor networks exchanging 
time stamps based on the time at their respective clocks. 
A survey of the popular approaches employed in practice 
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for timing synchronization is presented in [2] and [3|. The 
one-way message exchange mechanism involves a reference 
node broadcasting its timing information to other nodes in 
a network. The receiver nodes record the arrival of these 
messages with respect to their own clock. After several such 
time stamps have been exchanged, the nodes estimate their 
offsets based on these observations. A particular case of this 
approach is the flooding time synchronization protocol (FTSP) 
101 which uses regression to estimate the clock offset. On 
the other hand, through a two-way timing exchange process, 
adjacent nodes aim to achieve pairwise synchronization by 
communicating their timing information with each other. After 
a round of N messages, each nodes tries to estimate its 
own clock parameters. The timing-sync protocol for sensor 
networks (TPSNs) Q uses this strategy in two phases to 
synchronize clocks in a network. The level discovery phase 
involves a spanning tree based representation of a WSN 
while nodes attempt to synchronize with their immediate 
parents using a two-way message exchange process in the 
synchronization phase. In receiver-receiver synchronization, 
nodes collect time stamps sent from a common broadcasting 
node and utilize them to adjust their clocks. The reference 
broadcast synchronization (RBS) protocol [6| uses reference 
beacons sent from a master node to establish a common 
notion of time across a network. An alternative framework 
for network-wide distributed clock synchronization consists 
of recasting the problem of agreement on oscillation phases 
and/or frequencies as a consensus based recursive model in 
which only local message passing is required among nodes. 
By assuming a connected network, it is possible to design ef- 
ficient distributed algorithms by carefully choosing the update 
function. Under this framework, |9| proposed a Laplacian- 
based algorithm for establishing agreement on oscillation 
frequencies all over the network based on standard consensus. 
A combined agreement over both clock phases and frequencies 
has been studied in [10], by making use of state-of-the-art fast 
consensus techniques. Scalable synchronization algorithms for 
large sensor networks are developed in [11] and [12] inspired 
by mathematical biology models justifying synchrony in the 
biological agents. 

The clock synchronization problem in a WSN offers a 
natural statistical signal processing framework whereby, the 
clock parameters are to be estimated using timing information 
from various sensors l22l . A model based synchronization 
approach to arrest the clock drifts is explored in [8|. The 
impairments in message transmission arise from the various 
delays experienced by the messages as they travel through 
the transmission medium. Therefore, a crucial component of 
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efficient clock parameter estimation is accurate modeling of 
the network delay distributions. Several distributions have been 
proposed that aim to capture the random queuing delays in 
a network [7|. Some of these candidate distributions include 
exponential, Weibull, Gamma and log-normal distributions. 
Assuming an exponential delay distribution, several estimators 
were proposed in lfT3ll . It was argued that when the propagation 
delay d is unknown, the maximum likelihood (ML) estimator 
for the clock offset 9 is not unique. However, it was later 
shown in |14| that the ML estimator of 9 does exist uniquely 
for the case of unknown d. The performance of these estima- 
tors was compared with benchmark estimation bounds in 1 15 1. 
Considering an offset and skew model, Chaudhari et.al. pre- 
sented algorithms for the joint ML estimation of clock offset 
and skew in IT6l when the network delays are exponentially 
distributed. Clock offset and skew estimators were determined 
in ifTTl based on the assumption that the network delays 
arise from the contribution of several independent processes 
and as such, were modeled as Gaussian. The convergence 
of distributed consensus time synchronization algorithms is 
investigated in lfl8ll and fl9l . assuming a Gaussian delay 
between sensor nodes. More recently, the minimum variance 
unbiased estimator (MVUE) for the clock offset under an 
exponential delay model was proposed in ||20l . The timing 
synchronization problem for the offset-only case was also 
recast as an instance of convex optimization in |2D for Weibull 
distributed network delays. A recent contribution [23] has 
investigated the feasibility of determining the clock parameters 
by studying the fundamental limits on clock synchronization 
for wireline and wireless networks. 

In this work, considering the classic two-way message 
exchange mechanism, a unified framework for the clock offset 
estimation problem is presented when the likelihood function 
of the observations is Gaussian, exponential or log-normally 
distributed. A parameterized solution is proposed for the 
ML estimation of clock offset by recasting the likelihood 
maximization as an instance of convex optimization. In order 
to incorporate the effect of time variations in the clock 
offset between sensor nodes, a Bayesian inference approach 
is also studied based on a factor graph representation of the 
posterior density. The major contributions of this work can be 
summarized as follows. 

1) A unified framework for ML estimation of clock offset, 
based on convex optimization, is presented when the 
likelihood function of the observations is Gaussian, 
exponential and log-normally distributed. The proposed 
framework recovers the already known results for Gaus- 
sian and exponentially distributed likelihood functions 
and determines the ML estimate in case of log-normal 
distribution. Hence, the proposed convex optimization 
based approach represents a simpler alternative, and a 
more general derivation of ML estimator, which by- 
passes the graphical analysis used in lfl4l to maximize 
the likelihood function. 

2) In order to capture the time variations in clock offsets 
due to imperfect oscillators, a Bayesian framework is 
presented by considering the clock offset as a random 



Gauss-Markov process. Bayesian inference is performed 
using factor graphs and the max-product algorithm. The 
message passing strategy yields a closed form solution 
for Gaussian, exponential and log-normally distributed 
likelihood functions. This extends the current literature 
to cases where the clock offset may not be deterministic, 
but is in fact a random process. 
3) In order to evaluate the performance of the proposed 
estimators, classical as well as Bayesian bounds are 
derived for arbitrary exponential family distributed like- 
lihood functions, which is a wide class and contains 
almost all distributions of interest. While these results 
aid in comparing various estimators in this work, they 
can be useful in their own right in classical and Bayesian 
estimation theory. 

This paper is organized as follows. The system model is 
outlined in Section [II] The ML estimation of clock offset 



based on convex optimization is proposed in Section III The 
factor graph based inference algorithm for the synchronization 



problem in a Bayesian paradigm is detailed in Section IV and 
a closed form solution is obtained. Section [V] presents several 
theoretical lower bounds on the variance of an estimator eval- 
uated in the classical as well as Bayesian regime. Simulation 
studies are discussed in Section [Vl] which corroborate the 



earlier results. Finally, the paper is concluded in Section VII 
along with some directions for future research. 



II. System Model 

The process of pairwise synchronization between two nodes 
S and R is illustrated in Fig. [T] At the jth message exchange, 
node S sends the information about its current time through 
a message including time stamp Tj. Upon receipt of this 
message, Node R records the reception time Tj according to 
its own time scale. The two-way timing message exchange pro- 
cess is completed when node R replies with a synchronization 
packet containing time stamps Tj and Tj which is received 
at time Tj by node S with respect to its own clock. After N 
such messages have been exchanged between nodes S and R, 
node S is equipped with time stamps {Tj ,Tj ,Tj ,Tj}f =1 . 
The impairments in the signaling mechanism occur due to a 
fixed propagation delay, which accounts for the time required 
by the message to travel through the transmission medium, 
and a variable network delay, that arises due to queuing 
delay experienced by the messages during transmission and 
reception. By assuming that the respective clocks of nodes S 
and R are related by Ca(t) — 9 + Cs{t), the two-way timing 
message exchange model at the jth instant can be represented 
as 



rpl rril 

3 ~ 3 

rpA rp3 

3 ~ 3 



d + 9 + Xj 
d-9 + Y, 



(1) 



where d represents the propagation delay, assumed symmetric 
in both directions, and 9 is offset of the clock at node R 
relative to the clock at node S. X, and Yj are the indepen- 
dent and identically distributed variable network delays. By 
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d + e + x. 



3 



d-6 + Yi 



V = [Vt,. 



, Vat] is Gaussian or log-normally distributed 



is given below. 



Unconstrained Likelihood: 



/(U;Oocexp 



3=1 
N 



/(V; V) cx exp U £ ^ (Vj) - (iP) 



(5) 



(6) 



Fig. 1. A two-way timing message exchange mechanism 
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77- ^ T 2 



T 1 

3 
3 



the system in ([T} can be equivalently expressed as 

fl, 



By further defining 



IjAd- 



(2) 



(3) 



the model in |2) can be written as 

U j = ^ + X j 
V j = iP + Y j 

for j = 1,...,N. The goal is to determine precise estimates 
of £ and i/j using observations {Uj, Vj}^ =l . An estimate of 9 
can, in turn, be obtained using <|3j as follows 

= ^ ■ (4) 

Accurate modeling of the variable delays, X, and Yj, has 
been a topic of interest in recent years. Several distributions 
have been proposed that aim to capture the random effects 
caused by the queuing delays [7|. These distributions include 
exponential, gamma, log-normal and Weibull. In addition, the 
authors in [ 17 1 argued that Xj and Yj result from contributions 
of numerous independent random processes and can, therefore, 
be assumed to be Gaussian. The ML estimate of d and for 
the case of exponential distribution was determined in |14|. 
Recently, the minimum variance unbiased estimate (MVUE) 
of the clock offset under an exponentially distributed network 
delay was proposed in EDI . In this work, instead of working 
with a specific distribution, a general framework of the clock 
synchronization problem is proposed that yields a parameter- 
ized solution of the clock offset estimation problem in the 
classical as well Bayesian regime when the likelihood function 
of the observations, Uj and Vj, is Gaussian, exponential or 
log-normally distributed. 

In particular, the general notation used when the likelihood 



where r]^(Uj) and rj^(Vj) are sufficient statistics for estimating 
£ and ip, respectively. The log -partition functions </>f(.) and 
(f>tj,(.) serve as normalization factors so that /(U;£) and 
/(V; ip) are valid probability distributions. The likelihood 
function is called 'unconstrained' since its domain is inde- 
pendent of the parameters £ and ip. 

Similarly, the general notation used for an exponentially 
distributed likelihood function is given below. 

Constrained Likelihood: 

(N \ N 



(7) 



N \ N 

/(V; V) cx exp ( ^ £ r)^(Vj) - N<f>^) [] l(V 3 - V) 

J = l / 3=1 



(8) 



where the indicator function I(.) is defined as 

'l x > 



l(x) 



x < 



and the roles of r)^(Uj), r]^{Vj), (/>((■) and 4>^{.) are similar to 
(|5]l and |6]). The likelihood function is called constrained since 
its domain depends on the parameters £ and ip. It must be noted 
that the likelihood functions <|5j-<[8j are expressed in terms of 
general exponential family distributions. This approach helps 
to keep the exposition sufficiently general and also allows us 
to recover the known results for the ML estimation of clock 
offset for Gaussian and exponentially distributed likelihood 
functions |14] [15|, and determine the ML estimator of the 
clock offset in case of log-normally distributed likelihood 



function, as shown in Section III The proposed approach will 
also prove useful in investigating a unified novel framework for 
clock offset estimation in the Bayesian setting for Gaussian, 
exponential or log-normally distributed likelihood functions, 
as will be shown in Section [Tvl 

Some key ingredients of the proposed solution for the 
clock offset estimation problem, based on the properties of 
exponential family, can be summarized as follows ||25l . 

1) The mean and variance of the sufficient statistic r]^{Uj) 
are expressed as 



Efo C ([/;)] = 



function of the observations U 



Pi 



, Un] and 



4=Var[^)] 



(9) 



(10) 
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2) The moment generating function (MGF) of the statistic 
f]^(Uj) is given by 

M llt {h) = expfafc + h) -MO) ■ (ID 

3) The non-negativity of the variance cr^ in ( flO) implies 
that the log-partition function </>{(.) is convex. 

4) For Gaussian, exponential and log-normally distributed 
likelihood functions, the log-partition function <f)£ (£) can 
be expressed as a second degree polynomial given by 



(12) 



The coefficient in this approximation can be obtained 
using the variance of the statistic rj^(Uj), which is 
assumed known. Using (jTUJ, is given by 

< 

If the statistical moment in (JTOf is not available, the 
empirical moment can be substituted since it readily 
follows from the weak law of large numbers that 



N 

9/ 



E[^(t/)], iV->-oo 



Similar expressions can also be written for r]^(Vj), M v ^(h) 
and (f>^,(ip), respectively. Using the aforementioned properties 
of the exponential family, the ML as well as Bayesian esti- 
mates of £ and ip are to be determined utilizing the data set 
{Uj,Vj}^L v based on ([12}. 

III. Maximum Likelihood Estimation 

In this section, the ML estimates of 9 are obtained by 
recasting the likelihood maximization as an instance of convex 
optimization. This approach differs from the graphical argu- 
ments used to maximize the likelihood in |14|. The specific 
cases of unconstrained and constrained likelihood functions 
are considered separately. 

A. Unconstrained Likelihood 

Using |5]), |6} and ( |12) , the unconstrained likelihood func- 
tions are given by 



N 



/(U;0 cx exp UX>(^) 

( N a 2 
/(V;V)cxexp U]T>^)-A^V> 2 

The ML estimates of £ and ip can now be expressed as 

( N a 2 
£ml = arg max exp £2j»fc(I7 3 -) - N-^S? 

N _2 



(13) 



(14) 



V^ml = arg max exp ^y^V^jVj) - N-^ip 2 . (16) 



Theorem 1: The likelihood maximization problems ( fl"5j ) 
and ( [To} are strictly concave and the ML estimates are given 
by 



N < 



(17a) 



(17b) 



Hence, the ML estimator 9 ML for the clock offset can be 
written as 



7 ML — 



(18) 



Proof: The ML estimate of £ in ( fT5| can be equivalently 
determined by maximizing the exponent in the likelihood 
function. It can be easily verified that the exponent in ( |T5| , 
which is a quadratic function of £, is strictly concave |29|. 
A similar explanation applies to the ML estimate of ip. The 
ML estimates in ( fT7] i can be obtained by setting the first 
derivative of the exponent with respect to £ (resp. ip) to zero. 
The estimate 9ml in ( fl"8| > can be obtained by invoking the 
invariance principle [28 1. Hence, the proof readily follows. ■ 
1 ) Gaussian Distributed Likelihood Function: A particular 
application of Theorem 1 is the case when the likelihood 
functions /(£/?;£) and f(Vj]il>) have a Gaussian distribution 
i.e., f(U f ,0 ~ A^,af) and f(Vf^) ~ Af(ip,a 2 ) E). 
Hence, it follows that 



/(U;0 = 



1 



2a 2 



(2naf) 



exp 



(19) 



which can be rearranged as 



/(U;0«expU 



a 2 2a 2 



By comparing with ( fT3~| , we have 



and the ML estimate using ( |17a| > is given by 
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ML 



TV 



(20) 



(21) 



By a similar reasoning, the ML estimate for tp ( |17b| > is given 
by 



ML 



TV 



Using ( fT8| >, the ML estimate for the offset 9 can be expressed 



(15) as 



2N 



(22) 



The above estimate coincides exactly with the one reported in 

ma. 
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2) Log-Normally Distributed Likelihood Function: Another 
application of Theorem [T] is when the samples Uj and Vj are 
log-normally distributed. In this case, we have 



/(U;0 



=n^ iex p (- 

,JV 



2af 



oc exp £ 



cr| 2af 



(23) 



A comparison with ( fT3| l yields 

%(^) 



Following a similar line of reasoning, 

nr \ lo S V 3 2 



The ML estimator for 9 can obtained from ( fT~8] > using ( |17a[ ) 
and ( |17b| i, and is given by 



Ef=i (log^-log^) 



N 



(24) 



B. Constrained Likelihood 

Using Q, (|8]l and ( fT2| ), the constrained likelihood functions 
are given by 



JV 



Jl \ N 



/(U; Ocx exp - N^ 2 ]J l{Uj - £) 



3=1 



JV 



3=1 



,2 \ JV 



(25) 



/(V;^)cxexp ^^(yjJ-tf-^V ]J - V) 



7=1 



3=1 



(26) 



The resulting ML estimates of £ and -0 can be obtained as 



Iml = arg max exp ^^rj^Uj) - JV-|=£ 2 

^ 3=1 



such that Uj > £ 



N 



■0ml = arg max exp ^^(^) - N-^ip 2 



3=1 



such that V",' > 



(27) 



(28) 



Theorem 2: The likelihood maximization problems (J27j) 
and ( |28| are strictly concave and the ML estimates can be 
expressed as 



^ ML = min l v NaJ ( ^\ 

0ml = mm — 2 , V ( i) 



(29a) 



(29b) 



where Un\ and Vm denote the first order statistics of the 
samples Uj and Vj, respectively. The ML estimator 9ml for 
the clock offset is given by 



^ML 



7 ML 



(30) 



Proof: By using arguments similar to Theorem 1 and 
noting that the TAT constrains Uj > £ (resp. Vj > 0) are linear 
functions of £ (resp. ?/>), the proof of concavity readily follows. 
Also, the likelihood maximization problems ( |2T| and ( |28j ) can 
be equivalently expressed as 



A'* 



|ml = arg max exp €^Trid U j) " 



3 = 1 



such that Um > £ 



JV 



■0ml = arg max exp V>y)tfy.(Vj) - N^0 2 



3 = 1 

such that Vqj > ^ . 

The unconstrained maximizer £ of the objective function in 
p7| ) is given by ( |17a| ). If £ < f7( X ), then | ML = I- On the other 
hand, if Um < £, then the ML estimate is given by £ml = 
[/(i) using concavity of the objective function. Combining the 
two cases, the ML estimate in ( |29a| i is obtained. A similar 
explanation applies to the ML estimate 0ml in |29b) and 9ml 
follows from the invariance principle. ■ 
1 ) Exponentially Distributed Likelihood Function: For the 
case when the likelihood functions are exponentially dis- 
tributed, the density function of the samples Uj can be written 
asEl 

/(U;£)=Af exp f-A^^-tfj I(U (1) - £) (31) 

where AjT 1 is the mean of the delays Xj. The density function 
can be rearranged as 

/(U;£)c<exp(iVA e O . 

Comparing the above formulation with ( |25| ), 

Ve (Uj) = A 6 <r* = 0. (32) 

Using ( |29a| i, the ML estimate is given by 

Iml = U (1) . (33) 

Employing a similar reasoning, 

0>ML = V(i) . 

Using ( f30l >, the ML estimate of # is given by 

2 _ U{i) - V(i) 

"ml — ^ ' ^ -* 

which coincides exactly with the one reported in iTfll . where 
it is derived using graphical arguments. 

Remark 1: The ML estimation method outlined above dif- 
fers from the previous work 1 14 1 in that it is based on convex 
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optimization, while [14] maximized the likelihood graphically. 
Hence, this approach presents an alternative view of the ML 
estimation of the clock offset. It also allows us the determine 
the ML estimator of 9 when the likelihood function is log- 
normally distributed. In addition, Theorem [2] will also be 
useful in Section IV where estimation of 9 in the Bayesian 
regime is discussed . 



IV. A Factor Graph Approach 

The imperfections introduced by environmental conditions 
in the quartz oscillator in sensor nodes results in a time-varying 
clock offset between nodes in a WSN. To cater for such a 
temporal variation, in this section, a Bayesian approach to 
the clock synchronization problem is adopted by representing 
the a-posteriori density as a factor graph. The inference is 
performed on the resulting factor graph by message passing 
using max-product algorithm. To ensure completeness, a brief 
description of factor graphs and the max-product algorithm is 
provided below. 

A factor graph is a bipartite graph that represents a factor- 
ization of a global function as a product of local functions 
called factors, each factor being dependent on a subset of 
variables. Factor graphs are often used to produce a graphical 
model depicting the various inter-dependencies between a 
collection of interacting variables. Each factor is represented 
by a factor node and each variable has an edge or a half-edge. 
An edge connects a particular variable to a factor node if and 
only if it is an argument of the factor expressed by the factor 
node 1241 . 

Inference can be performed by passing messages 
(sometimes called beliefs) along the edges of a factor 
graph. In particular, max-product algorithm is used to 
compute the messages exchanged between variables and 
factor nodes. These messages can be summarized as follows 

variable to factor node : 



l x ^f (x) = Y\_ m h-+x (x) 
h£n{x)\f 



(35) 



factor node to variable 



max 



f(z) n 



(36) 



zGn(/)\{x} 



where Z = n (/) is the set of arguments of the local function 
/. The marginal distributions associated with each variable can 
be obtained by the product of all incoming messages on the 
variable. 

In order to sufficiently capture the temporal variations, the 
parameters £ and ip are assumed to evolve through a Gauss- 
Markov process given by 



6c — £fe-i - 
ip k = ipk-i 



- v k 



for k = 1, 



,N 



where Wk and Vk are i.i.d such that Wk,Vk ~ A/"(0, er 2 ). The 
posterior pdf can be expressed as 



m,ii>\u,v)<xf(s,ij>)f(u,v\e,ii>) 



N 



N 



/&>) n /(fki&-i)/(Vto) n fWkWk-i) 



*:=i 



fc=i 



N 



Hf(Uk\z k )f(v k m 



(37) 



k=l 



where uniform priors /(£o) and f(ipo) are assumed. Define 

M{A-i,a 2 ), fk = f(U k \t k ), h k = f(V k \i>k), where the 
likelihood functions are given by 



/(C/ fc |£ fc ) cx exp l£ k T,t(U k ) - 



f<y k \ik) cx exp \^ k Vi,(V k ) - -j±i> 2 k j , (38) 

based on (12) . The resulting factor graph representation of the 
posterior pdf is shown in Fig. [2] 

Remark 2: A few important observations of this represen- 
tation can be summarized below. 

• Notice that the substitution in ^ renders a cycle-free na- 
ture to the factor graph. Therefore, inference by message 
passing on such a factor graph is indeed optimal ll24ll . 

• The two factor graphs shown in Fig. 2 have a similar 
structure and hence, message computations will only be 
shown for the estimate £at. Clearly, similar expressions 
will apply to ipN- 

In addition, only the case of constrained likelihood will be 
considered, since the case of an unconstrained likelihood is 
subsumed, as will be shown shortly. The clock offset estimator 
9n can be obtained from £jv and ip^ using Q. 

By defining = — ^ and ,% fe = rj^(U k ), the 

constrained likelihood function for the samples U k can be 
written as 



f k (x exp (a S}k £,l + /8£,*£fc) l{U k - f fc ) 



(39) 



The message passing strategy starts by sending a message 
from the factor node fa to the variable £jy The variable 
£jv relays this message to the factor node S^_ v The factor 
node computes the product of this message with the factor 
and sends the resulting message to the variable {jjy-i 
after 'summarizing' over the variable £jv. In the max-product 
algorithm, a 'max' function is used as a summary propagation 
operator (cf. ([361)). These messages are computed as 

m /jv^?« = In 



cx max S 



N 

N-l 



max : 
5« V27r(T- 



: exp 



-(6v — €n-iY 

2fi 2 

exp (a ( , N £,ff + /Sf.jv^jv) K~U N - £ N ) 
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"N-l 
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i>i 
















v° 
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v k-l 



V N-1 



Fig. 2. Factor graph representation of the posterior density l |37) 



which can be rearranged as 



m S N _ >5jv _ 1 oc max exp( J 4 ?>A r^ + B 5iA r£w-i + 



where 



. A 1 A 1 



(41) 



Let (jv be the unconstrained maximizer of the exponent in the 
objective function above, i.e., 

£at = argmax (A itN $,% + -Bj./vCat-i + Q,/v6v£/v-i + 



This implies that 

6v 



2A 



(42) 



Following a line of reasoning similar to Theorem 2, it follows 
that 

£ N = min (fjv, J/jv) ■ 

However, £jv depends on £/v-i> which is undetermined at this 
stage. Hence, we need to further traverse the chain backwards. 
Assuming that £jy < Ujy, £jy from ( |42] i can be plugged back 
in (Hob which after some simplification yields 



(43) 



m ^_!^-i « ex P I ( ~ t¥ L \ 
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2A 



N-l 



The message passed from the variable £/v-i to the factor node 
$n-2 * s m e product of the message |43) and the message 
received from the factor node /jv-i, i-e-, 

S Jf-l^W-l /JV-1-+W-1 



Upon receipt of this message, the factor node (^I^ delivers 
a product of this message and the factor S^Z^ t0 the variable 
node £jv_2 after summarizing over £jv-i- This message can 
be expressed as 



m X N-i t , oc 



max 

£n-i<Un-i 



sN—l 
J N-2 



1 



=max —j= 
V27r(7 

• exp 



: exp 



AT-l 



5, 



r? 2 



4A 



2ct 2 



t 2 C^nD^n c 

?JV-1 ^1 



2A 



• exp (a^,jv_i^_ 1 + /%jv-i6v-i) HUn-i ~ 6v-i) ■ 

After some algebraic steps, the message above can be com- 
pactly represented as 



™ 5 « i _ cx max exp(A CN ^ N _ 1 + 

B£,N-i£n-2 + C?,Af-lCjV-l^JV-2 + ^S.JV-lCiV-l) (44) 

where 



1 

2a 



2 + a£,jv-i + B^.n 



r 2 
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Be, N-l 



2a 



,N-1 



1 r A 1 

2^,w ' 



/v-i 



Proceeding as before, the unconstrained maximizer £at_i of 
the objective function above is given by 



AT-l 



and the solution to the maximization problem (j44j» is expressed 

£at_i = min (fjv-i, t/jv-i) • 



as 



Again, £jv_i depends on £jv_2 and therefore, the solution 
demands another traversal backwards on the factor graph 
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representation in Fig. 2. By plugging £jv-i back in ( |44| i, it 
follows that 



m<-N-i oc 



exp < Bt tN -! 



N-l \ ^2 C^.JV-l-Df.jV-l 



4A 



£,AT-1 



(45) 

which has a form similar to fOj ). It is clear that one can keep 
traversing back in the graph yielding messages similar to ( |4"3"j ) 
and (|45]>. In general, for i = 1, . . . , N — 1 we can write 



(46) 



i r 2 

i A 1 i id °£,iV-i+X 

2<H 4Ai,Ar_, i+1 



A 



6,N-i 



1 

2a2 



A 



n A o C^JV-i+i-D^JV-i+l 



2A 



and 



t Q.JV-ifiV-i-l + -Df.AT-i 

4jv-j ^ (4/) 

£,N~i = mm(£ N _i,U N -i) ■ (48) 
Using ((47]) and (|48]> with i = iV - 1, it follows that 
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ll = min (£l>^l) 



(49) 
(50) 



Similarly, by observing the form of ( |4"3"j ) and ( |45| , it follows 
that 



a 2 



(51) 

The estimate £ can be obtained by maximizing the received 
message in < f5T~| >. It can be noticed from the structure of the 
factor graph that this maximization is inherently unconstrained 
i.e., 

lo = lo = max m s i_+* 
so 



(52) 



The estimate in ( |52| ) can now be substituted in ( |49] i to yield 
£i, which can then be used to solve for £i in |50|. Clearly, 
this chain of calculations can be continued using recursions 
<g7) and (g|). 
Define 

a Q,fex + -D{,fc 



2A 



(53) 



A key property of the function <?£,&(.), which proves useful 
in the quest for a closed form solution, can be summarized in 
the following lemma. 

Lemma 1: For real numbers a and b, the function 
defined in d53l satisfies 



g^ k (min(a,6)) = min (g e>fc (a), fl£,fc(6)) 



Proof: The constants A^fc, C^.fc and are defined in 
< |4TT > and |46|. The proof follows by noting that > 

which implies that <7£,fe(.) is a monotonically increasing func- 
tion. ■ 
With the notation <7f,fc(.)> the following chain of equalities can 
be conveniently written as 

li = S4U (lo) 

ii = min (Ui,9i,i (lo)) 

6 = 5«,2 (|i) 

| 2 = min (V 2 ,5 5 ,2 (|i)) 

where 

.95,2 (li) = .%,2 (min (lJi,g^i (lo))) 

= min ( 5c , 2 (E7i) ,#, a (54,1 (lo))) (54) 

where ( |54| i follows from Lemma [T] The estimate | 2 can be 
expressed as 

| 2 = min (u 2 ,mm (g ii2 (Ui),g^ 2 (95,1 (lo)))) 
= min (u 2 ,g ij2 (^1) ,3£,2 (flti (lo))) ■ 
By following the same procedure, one can write 



£3 = min (^C/ 3 ,5e,3 (U 2 ) ,3£,3 (Sf,2 (^1)) , 
3C3 (54,2 (54,1 (lo)))) • 



For m > j, define 



A 



G&(0=$f,m(#,m-l...«fej(.)) • (53 > 

The estimate £3 can, therefore, be compactly represented as 

Is = min (lf 3 , Gf i3 (C/ 2 ) , Gf, 2 (C/i) , G|,a (| ) ) . 

Hence, one can keep estimating at each stage using this 
strategy. Note that the estimator only depends on functions of 
data and can be readily evaluated. 

In order to derive analogous expressions for ip, a similar line 
of reasoning should be followed. In particular, the constants 
A CiA r_i, B 6tN -i, C^N-i and D £iN _i for i = 0, . . . , N - 1, 
can be obtained straightforwardly from ( |4Tj ) and |46} by 
substituting a^-i an d /?£ jv-< with a 4>,N-i an d Pib,N-ii 
respectively. Using these constants, ipo, g^^ a nd ■ can 
be defined analogously to ( |52") >, ( |53"] > and ( f55| . 

Generalizing this framework, the closed form expression for 
the clock offset estimate 6m is given by the following theorem. 

Theorem 3: The state estimates £jy and ip^ for the posterior 
pdf in ( (37} can be expressed as 

|jv = min (t/jv, G£ N (f/jv-i) , • ■ ■ , G^ 2 (f/i) , (| ) ) 
^ = min (Ujy, G^ iW (^v_i) , . . . , G^ 2 (Vx) , Gf a (^ ) ) 
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and the factor graph based clock offset estimate (FGE) On is 
given by 



6v - ^ 



N 



(56) 



Proof: The proof follows from the discussion above and 
using Q. ■ 

Remark 3: The closed form expressions in Theorem [3] 
enable the estimation of the clock offset, when it may be time- 
varying and the likelihood functions, f(Uk\£,k) and f(Vk\il>k), 
have a Gaussian, exponential or log-normal distribution. 



A. Gaussian Distributed Likelihood Function 

A particular case of the Bayesian framework described 
above occurs when likelihood functions f(Uk\£,k) and 
f{Vk\ipk) have a Gaussian distribution, i.e., f(Uk\£k) ~ 
JV(&,of >fc ) and /(Wfc) ~-Wk><*). Le - 



f{Uk\Zk) 



exp ■ 



UkUk 



2cr 



JL 



(57) 



The aforementioned Gaussian distribution constitutes an un- 
constrained likelihood function, i.e., the domain of the pdf is 
independent of the unknown parameter It is clear from the 
message passing approach that at each stage k of the factor 
graph, the unconstrained maximizer ^ is the actual solution 
to the likelihood maximization problem 



max exp (A^l + Sf,fcCfc-i + C i,ktk€k-i + D^ k ) 



i.e., £k = £,k Vfc = 1,...,N. Hence, the unconstrained 
likelihood maximization problem is subsumed in the message 
passing framework for constrained likelihood maximization. 
It follows from Theorem [3] that £n f° r Gaussian distributed 
observations Uk in d57) is given by 



£n - G^i (£0) 



where £0 and G^(.) are defined in ([52]) and ( |55] l, respectively. 
Evaluating £jy requires to determine the constants in ( |4Tj ) and 
(|46]i. By comparing (j57j with |39|, we have 



2 °h ' 



(58) 



Using these values for a^k and /3f,fe> < FT) and |46) can be 
written as 

1 1 B - 1 



2a 2 



1 



2a 2 



£7 



JY 



(59) 



1 

2^2 

1 

2^' 



£,JV-i+l 



r 2 



4A 



•£,iV-t+l 



(T 



2A 



£,JV— i+l 



, JV— 1, Using similar arguments, it can be shown 



for i = 1, 
that the estimate -i/V is given by 

= ($>) ■ 

It follows from Q that the FGE, On, can be expressed as 



' N = 



(60) 



The behavior of On can be further investigated for the case 
when the noise variance a 2 in the Gauss-Markov model goes 
to zero. Consider 

9t,N{£) = tt& 

where the constants A^n< ^£,at> C^jv and n are given by 
( |59l >. After some algebraic steps, we have 

As a 2 — > 0, 5f.jv(£) £■ Similarly, it can be shown that 
9£,N-i(0 — > 6 as a 2 — >• 0. Hence, it follows that in the low 
system noise regime, as cr 2 — >• 

Similarly, it can be shown that 



Therefore 



AA^^B^^i — C 2 ! 
£0 ~ $0 



which can be proven to be equal to the ML estimator ( |22j i. 

B. Log-Normally Distributed Likelihood Function 

The log-normally distributed likelihood function in the 
Bayesian regime can be expressed as 

\2\ 



f(u k \£ k ) = 



1 



: exp 



(iogt/ fc -ar 



UkCr^W^n 

'z k io g (u k ) e k 



oc exp 



(61) 



10 



By comparing ( |6T| i and p9] l, we have 

1 „ logC/ fc 

£,k £.* 
Clearly, the only difference here with the Gaussian distribution 
is a redefinition of (3^ The expression of £jv in this case is 
again 

where G^ 1 (.) and £ are given by ( |55] l and ( [52^ , respectively. 
The recursively evaluated constants in ( |4T| > and |46} can be 
written as 



logf/Af 



2a 2 2<t| !JV _. 



Bf.N-i — —-. 



1 



1 



log?7Ar_i Cf^-i+l-DgjJV-i+l 
A , - ., - 



0" 



£,N-i 



2^5,jv-i+i 



for i = 1, . . . , N — 1. Similar arguments apply to ^i. Hence, 
the FGE 9jf can be expressed as 



7/V = 



(lo)-G^ (^ ) 



(62) 



Again, as the Gauss-Markov system noise er 2 — > 0, the above 
estimator approaches its ML counterpart (|24]i. 



C. Exponential Distribution 

Theorem [3] can also be used to derive a Bayesian estimator 
£jv for the exponentially distributed likelihood case considered 
in lfl4ll . In this case, we have 

f(Uh\tk) = A e cxp (-X^Uk - &)) I(C/ fe - &) 

ex exp(A £ a)I(C/fc - &) (63) 
where A^ 1 is the mean network delay of Xk- A comparison 



of ( |63| l with p9f reveals that 

a ?,fc = 0, fi^u = A ? . 

For these values of a^fc and /3{,fe, the constants A^k, B^^, 
C^.fc and ZJ^ fc are given by 

for all k = 1, . . . , AT. Using Theorem [3] we have 



G^ N (U N -i) - - 



2A 6tN 

2 



= U N -i + X^a- 
Similarly it can be shown that 

G? n _ 1 (Un-2) — Un-2 + 2A^(T 2 



and so on. The estimator £o at the last step can be evaluated 

Co = - G|, = + °° ' 



This implies that 



= +oo 



(64) 



Using (|64]i and Theorem [3] it readily follows that 

|jv = mm(U N , U N -i + A ? cr 2 , U N -2 + 2A 4 ct 2 , 

...,C/i + (JV- l)A c a 2 ) . (65) 

Similarly, for the pdf OH 

/( W fc ) = X 4 , exp (-A^(U fe - Vfc)) I(Vfc - 
oc exp(A^,-0fc)I(Vfc - r/'fc) , 

the estimate t/>at is given by 

^iv = min(Vjv, Vn-i + X^o 2 ,V N - 2 + 2A^ct 2 , 

. . . , Vx + {N - l)A^a 2 ) (66) 

and the estimate On can be obtained using ( |56] >, ( |65] > and ((66]) 

as 

0jv = ^ rnin(f/;v, L^-i + A ? <r 2 , U N -2 + 2A 5 cr 2 , 



1 



■ ■ ■ ,Ui + (N — l)Ag(T )— 
min(VAr, Fjv_! + A^ct 2 , V/v_ 2 + 2A^cr 2 , 



...,F! + (iV-l)A^<T 2 ) . (67) 
As the Gauss-Markov system noise a 2 0, ( |67| i yields 

- - min(C/jv,...,t/i) - min (V^, . . . , Vi) 
Wjv — > "ml — ^ ' 

which is the ML estimator given by ( f34j >. 

V. Classical and Bayesian Bounds 

To evaluate the performance of the estimators derived in the 
preceding sections, classical as well as Bayesian lower bounds 
on the variance of the estimators are discussed. The placement 
of a lower bound allows one to compare estimators by plotting 
their performance against the bound. It must be emphasized 
here that the results in this section assume no specific form of 
the log-partition function and are therefore, valid for arbitrary 
distributions from the exponential family, which is a wide class 
and contains almost all distributions of interest. Hence, these 
results are fairly general and can be useful in their own right 
in classical as well as Bayesian parameter estimation theory, 
and at the same time will be used as a stepping stone towards 
comparing the estimators developed thus far. 

The likelihood function of the data is considered an 
arbitrary member of the exponential family of distributions. 
In addition, depending on whether the domain of the 
likelihood depends on the parameter to be estimated, both 
cases of unconstrained as well as constrained likelihood 
functions are discussed to maintain full generality. The 
general expressions for the unconstrained and constrained 
likelihood functions for observations Z = [Z±, . . . , Zn] T are 
given by 
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Unconstrained Likelihood: 



N 



f(Z i p)<xew\pjyi(Z j )-N<t>(p) | i68s 



Constrained Likelihood: 



N 



N 



/(Z; p) ex exp pY.rii.Zj) - N<j>(p) J] l(Zj - p) (69) 

\ 3=1 / 3=1 

where p is the scalar parameter to be estimated. The goal 
is to derive lower bounds on the variance of estimators of 
p. For the case of classical estimation, the Cramer-Rao and 
Chapman-Robbins bounds are considered, while the Bayesian 
Cramer-Rao bound and a Bayesian version of the Chapman- 
Robbins bound are derived for the Bayesian paradigm. 

A. Cramer-Rao Lower Bound 

The Cramer-Rao lower bound (CRB) is a lower bound on 
the variance of an unbiased estimator of a deterministic pa- 
rameter. It is useful primarily because it is relatively simple to 
compute. However, it relies on certain 'regularity conditions' 
which are not satisfied by constrained likelihood functions 
when the domain of the likelihood depends on the unknown 
parameter (cf. ((69])). Hence, CRB is determined for the case 
of unconstrained likelihood functions only. 

In particular, CRB states that the variance of an unbiased 
estimator of p is lower bounded by 

Var(p) > ra9 lL..M ■ ( ? °) 



E 



d 2 ln/(Z;p) 
dp 2 



Theorem 4: The CRB for p in the unconstrained likelihood 
function in ( |6"8j ) is given by 

1 



Var(p) > 



Na 2 



(71) 



where 



2= d 2 <P(p) 

7 " dp 2 



Proof: The Fisher information for the likelihood function 
is given by 



I(p) £ E 



dp 2 



= -N^ = -No* 
dp 2 " 

and the proof readily follows. ■ 

B. Chapman-Robbins Bound 

Chapman-Robbins bound (CHRB), proposed in [26 1, sets a 
lower bound on the variance of an estimator of a deterministic 
parameter. The CHRB does not make any assumptions on the 
differentiability of the likelihood function and the regularity 
conditions that often constrain the use of CRB, and is substan- 
tially tighter than the CRB in many situations. Hence, CHRB 



is employed to determine a lower bound on the variance of an 
unbiased estimator of p for constrained likelihood functions. 
In general for a parameter p, the CHRB is given by 

'f(Z;p + hy 2 " 



Var(p) > 




/(Z;p) 



- l 



(72) 



which can be evaluated as shown below. 

Theorem 5: The CHRB for the parameter p given the 
likelihood function |69) can be expressed as 

-l 



Var(p) > 



{(M n (h))- 2N -C N (h)-l\ 



inf 

h 



h? 



where M n (h) is the MGF of the statistic rj(Zj) and 



CO) = E [exp (2hrj(Z j ))I(Z j -p-h)] 



(73) 



(74) 



with the expectation taken with respect to any Zj. 

Proof: The details of the proof are relegated to Appendix 
\M M 

C. Bayesian Cramer-Rao Lower Bound 

The Bayesian Cramer-Rao bound (BCRB) is a lower bound 
on the variance of an unbiased estimator when the parameter 
assumes a prior density. It requires the same regularity condi- 
tions to be satisfied as its classical counterpart. 

For an estimator p^ of p^, the BCRB states that the variance 
of the estimator is bounded below by the lower-right sub- 
matrix of the inverse of the Bayesian information matrix, 
Jcn(k) M, i.e., 

Var (p k ) > J c ^(k) = [J^(k)] kk 

where the Bayesian information matrix is given by 

~dlogf(Z k ,p k ) 01og/(Z*,p fc ) 



(75) 



[J C R(fc)]ij =E 



-E 



dpi 

d 2 logf(Z k , Pk ) 
dpidpj 



where the expectation is taken with respect to the joint pdf 
and 



Z k = [Z u . 



■ , z k 



Pk = [po,pl,---,Pk] 

f(Z k \p k ) ex exp (r)(Z k )p k - <p k {Pk)) 



(76) 



It is assumed that the parameter p k evolves through a Gauss- 
Markov model given by 

f(p k \ Pk -i) - ^cx P (- {pk - p r )2 ) (77) 



V2~T 



2a 2 



A recursive formula to evaluate the Bayesian sub-matrix, 
derived in [27 1, is given by 



Jcn(k + 1) = - E, 
+ E, 



(2) 
CR 

(3A) 
CR 



(k) (j CR (fc) 

(k)+E^\k) 



(78) 
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where 



and 



4 3 R B) (fc)^E 



•^-2 iog/(/? fe+1 |p fc ; 



y 2 



dp k dpk + i 
J 2 

d Pl+i 
d 2 



iog/(pfc+i|p fc ; 



j log/(pfc+i|pfe) 



log/(^fc+l|pfe+l 



and the expectation is again with respect to the joint pdf. 

Theorem 6: For the Bayesian framework in (|76j> and f77) , 
the recursive Bayesian information matrix in (|78j) is given by 



J C R(k + l) = (<J 2 + J c l(k)) 1 +o 



(79) 



with J C r(0) = 0. 

Proof: For the density functions, f(p k \p k ~i) an d 
/(■Z^felPfe) m f76), it can t> e verified that 



and 



(35) 
CR 



(fc) 



4to 



d 2 (j>k{pk+i) 
d Pl+i 



1 



(3A) 
CR 



(fc) = 



1 



d 2 <j) k (pk+i) 



d P 2 +l 

The proof follows by plugging these quantities in d78 



D. Bayesian Chapman-Robbins Bound 

A Bayesian version of the Chapman-Robbins bound 
(BCHRB) can be used to provide a lower bound on the 
variance of an estimator of p k when there are no regularity 
assumptions on the likelihood. In fact, unlike the BCRB, the 
BCHRB can be evaluated for constrained likelihood functions 
where the domain of the likelihood is dependent on the 
unknown parameter. 

BCHRB states that the variance of an estimator p k of p k 
is lower bounded as 

Var(p fc ) - [T fc (h fc ) - If 1 h fc h^ y 

with >z in the positive semi-definite sense, where 



T k (h k ) = E 



f(Z k ,p k + h k ) 



and h fc = [0, h lt . . . , h k ] T . 

Theorem 7: The BCHRB for the parameter p k can be 
expressed as 



Var(p fc ) > 



1 



where 



Jcn,k 
T k (h k ) - 1 



JcH,k — inf 2 



T k (h k )=\l[M- 2 (h J )M v (2h ] 



vj'=i 



exp 



1 2 



(80) 



Proof: See Appendix |B| for details. ■ 

Remark 4: Since the performance bounds are derived for 
arbitrary exponential family distributions, they can also prove 
useful in a broad sense in classical as well as Bayesian 
parameter estimation. 



E. Relation to Clock Offset Estimation 

The performance bounds derived in the preceding subsec- 
tions can be used to lower bound the mean square error 
(MSE) of the clock offset 8. Notice the similarity between 
the unconstrained and constrained likelihood functions for £ 
and %/} in |5|-([8]) and the general exponential family likelihood 
function considered in ( |68j ) and (69) . Therefore, the bounds 
derived above for p are also applicable to the parameter £ 
(and also ip). The MSE of the clock offset 9 can, in turn, be 
lower bounded using the bounds on £ and ip. 

Using Q, the following result is immediate. 

Proposition 1: The MSE of any estimator of 9 can be 
expressed as 



MSE 



(«) 



Var U 



Var 



<h 



where and are the biases of the estimators £ and ip, 
respectively. 

The explicit lower bounds on the MSE of any estimator 
9 for classical as well as Bayesian framework in case of 
Gaussian and exponentially distributed likelihood functions 
can be evaluated as shown below. 

1 ) Gaussian Distribution - CRB: If the likelihood function 
for £ is Gaussian distributed ( fT9| ), then using ( pO) and f7T) , 
it is straightforward to see that the CRB for any unbiased 
estimator £ is given by 



Var 



> 



N 



and a similar expression is applicable to ip as well. Using 
Proposition [T] it can be concluded that 



MSE 19) > 



4N 



(81) 



As a remark, it is evident in this case that 9 ML ( |22| ) is efficient 
in the sense that its MSE achieves ([81") with equality (cf. 



Appendix C-Ai 
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2) Exponential Distribution - CHRB: If the likelihood for 
£ is exponentially distributed 031) , using ( fTTj ) and ( |32j >, it can 
be easily verified that 

M V((u) (h) = 1 

and ( |74] > becomes 

C(/i) = exp (A 5 /i) , 
so that the statement of the CHRB ( f73] l can be rewritten as 



Var 



> 



inf 

h 



exp (A ? WV) - 1 



h 2 



0.6476 
ApV 2 " 



and similarly for tp. Using Proposition [I] it follows that 

(</<)) + 4 



MSE ( 9 ) = - { Var ( £ 



> 



0.162 

iV 2 




Var 



1 



(82) 



3) Gaussian Distribution - BCRB: In the Bayesian regime, 
if the likelihood function for £ is Gaussian distributed ( [57) 1, 
by using ( |58] > and ( |79| ), it can be seen that 



Jcr,c (A + 1) = (a 2 + (k) 



with Jcr,{ (0) = 0. A similar line of reasoning can be followed 
to derive an analogous recursion for Jcr,^. (k). The MSE of 
9 can be now be lower bounded as 



(83) 



4) Exponential Distribution - BCHRB: If the likelihood for 
£fc is exponentially distributed (|63"), (f80]> turns out to be 



A" 



Tfc(hfc) = exp [ A{ ^2 h o ) ex P 



1 fe 2 



In fact, we just have to notice that tfi^ is a constant 
function over and rj^(Uj) = Aj, so that ( f85] > becomes 



E [exp (2hjT]^(Uj))] = exp (X^hj 



therefore S'(hfc) = exp ( Af ^jLi h 



VI. Simulation Results 

This section aims to corroborate the theoretical results in 
preceding sections by conducting simulation studies in various 
scenarios. The performance of both the classical as well as 
Bayesian estimators is to be investigated. The measure of 
fidelity used to rate this performance is the MSE of the esti- 
mators for 9 and 9^. The parameter choice is = = 0.1 
for both Gaussian and log-normally distributed likelihoods, 
while A;; = A^, = 10 for exponentially distributed likelihood 
functions. 



X--X-X--X-X--X-X-X-X-X-X-X-X-X-X-X- 

* x X X 



x xx x x 



X X X X X X -0-CHRB 



-X— MLE 

- X - MLE Gaussian Based 
X MLE Exp Distr Based* 
- CRB 




Fig. 3. MSE and bounds for estimating 8 by using the MLE with log-normal 
likelihood. 



A. Log-normal Distribution 

No solution is reported thus far in literature in case the like- 
lihood is log-normally distributed. The proposed estimators, 
#ml ([24]) and 9n ( |62| ), can be used to determine an estimate 
for the clock offset in classical and Bayesian framework, 
respectively, as shown below. 

1) Classical Estimation Framework: The existing ap- 
proaches in literature only consider the Gaussian and the 
exponential cases, therefore ( [24] ) is a new result in the state- 
of-the-art about clock offset estimation. Fig. |3] shows a com- 
parison between the proposed ML estimator (MLE) p4] ) in 
case of a log-normally distributed likelihood ( |23] l with MLEs 
which (wrongly) assume that the likelihood is Gaussian and 
exponentially distributed, respectively. The plot shows that the 
latter approaches are not robust with respect to the likelihood 
distribution, and their performance is extremely poor if their 
assumptions do not hold. In addition, Fig. [3] also shows that 
the proposed MLE ( |24] i is efficient since it attains the CRB 
(as well as the CHRB). 

2) Bayesian Estimation Framework: Fig. |4] plots the MSE 
performance of the FGE ( |62| ) as well as the BCRB and 
BCHRB when the likelihoods are log-normally distributed 
( |6T] ), and a = 10~ 4 . Firstly, it can be seen that the MSE 
of the proposed FGE coincides with the estimation bounds. 
Secondly, as in the classical estimation case, if we were to 
(wrongly) assume a Gaussian or exponential distribution for 
the likelihoods ( [38] ), the resulting FGEs would perform poorly, 
a fact is evident in Fig.|4]by observing the unboundedness and 
unpredictability of the dashed curve (Gaussian assumption for 
the likelihoods) and the dotted curve (likelihoods assumed ex- 
ponentially distributed). This clearly establishes that the FGE 
([62]), obtained assuming log-normally distributed likelihoods, 
allows a strong performance improvement with respect to 
existing estimators if the likelihood functions ( |3~8] l are actually 
log-normally distributed. 
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Fig. 4. MSE and bounds for estimating 9^ by using FGE with log-normal Fig. 5. MSE and bounds for estimating 8 by using the MLE with Gaussian 



likelihood. 



and exponentially distributed likelihood. 



B. Estimator Performance vs Estimation Bounds 

It will also be useful to asses the performance of the 
MLEs in Gaussian and exponential cases derived in Section 
[Illjagainst the various benchmark estimation bounds derived in 
Section [V] Similarly, the FGEs for Gaussian and exponential 
distributions, proposed in Section IV can also be compared 
with the Bayesian bounds to study their MSE performance. 

1) Classical Estimation Framework: Fig. [5] shows the per- 
formance comparison between the MSE of the MLEs (|22j» 
and ([34} for Gaussian and exponentially distributed likelihood 
functions against the CRB and the CHRB, respectively. Firstly, 
it is evident that in the case of Gaussian distribution, the 
CRB and the CHRB coincide. Moreover, the MSE of §ml 
also coincides with the aforementioned bounds. On the other 
hand, for an exponentially distributed likelihood function, due 
to its lack of regularity, the CRB cannot be derived, thus only 
the CHRB is shown. It can be observed that the MSE of Oml 
is fairly close to CHRB, even though it does not coincide 
with it. From Fig. [5] the MSE of the MLEs for the Gaussian 
and exponential distribution case can be also compared. In 
order to ensure a fair comparison, parameters are chosen in 
a way to have the same variance of the observations for 
both distributions. From the MSE curves, one can infer that 
the MSE in case of an exponentially distributed likelihood 
is lower than the one for a Gaussian distribution as the 
number of observations N increases. This behavior is expected 
since it can be verified by the MSE expressions ([86} and 
([87}, reported in Appendix [C] that in case of a Gaussian 
distribution, the MSE decreases proportionally to 1/N, while 
in the exponential distribution case it decreases proportionally 
to 1/N 2 . 

2) Bayesian Estimation Framework: In Fig. [6j the MSE 
performance of the FGEs 6^ ( |60} and |67} is compared 
with BCRB and BCHRB for a = 1(T 4 . As in the classical 
estimation scenario, it is evident that for Gaussian distributed 
likelihoods, the MSE using ( |60} for 6n coincides with the 
reported bounds. The MSE of the FGE derived assuming 



— X- — FGE - Gaussian 
— B — BCRB - Gaussian 
BCHRB - Gaussian 

- X - FGE - Exp Distr 

- - BCHRB - Exp Distr 




Fig. 6. MSE and bounds for estimating t9jv by using FGE with Gaussian 
and exponentially distributed likelihood. 



exponentially distributed likelihoods ([67} is plotted against the 
BCHRB as well in Fig. [6] It is clear that the MSE is quite 
close to BCHRB, although not coinciding with it, as exactly 
was the case in the classical estimation framework. 

C. Comparing Classical and Bayesian frameworks 

The estimators proposed in the classical and the Bayesian 
framework can also be compared with each other based on 
their MSE performance as the system noise decreases. The 
aim here is to show that the latter approaches the former as 
a -> 0. 

Fig. [7] depicts the MSE for the cases of Gaussian, expo- 
nential and log-normal distribution for the likelihoods with 
N = 25. In the plot, the horizontal lines represent the MSEs 
in the classical framework, obtained with the MLEs, as shown 
in ( [86} and ( [87} in Appendix [C] It can be observed that, for all 
the three considered distributions, the MSE obtained by using 
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O Gaussian 
□ Exp- Distr 
X LogNormal 



B □ H B □■ — ■€> H B = - -U B B 



Fig. 7. MSE in estimating of 0jy vs a - Solid lines: FGE; dashed lines: MLE 
(classical framework). 



the FGE for estimating 8 approaches the MSE of the MLEs 
as cr < 10~ 3 . 



VII. Conclusions and Future Work 

The clock synchronization problem in sensor networks has 
received keen attention recently owing to its central role in crit- 
ical network operations as duty cycling, data fusion, and node 
localization. Based on a two-way timing message exchange 
scenario, this work proposes a unified framework for the 
clock offset estimation problem when the likelihood function 
of the observation time stamps is Gaussian, exponential and 
log-normally distributed. A convex optimization based ML 
estimation approach is presented for clock offsets. The results 
known thus far for Gaussian and exponentially distributed 
network delays are subsumed in the general approach while the 
ML estimator is derived when the likelihood function is log- 
normally distributed. In order to study the case of a possibly 
time-varying clock offset, a Bayesian approach is also studied 
using factor graphs. The novel message passing strategy results 
in a closed form solution of the time-varying clock offset 
estimation problem. In order to compare various estimators, 
several lower bounds on the variance of an estimator have 
been derived in the classical as well as the Bayesian regime 
for likelihood functions which are arbitrary members of the 
exponential family, a wide class containing several distribu- 
tions of interest. The theoretical findings are corroborated by 
simulation studies conducted in various scenarios. 

In future, it will be useful to incorporate the effect of 
clock skew in the clock offset estimation model. This can 
result in further reduction of the re-synchronization periods. 
In addition, the results about pairwise synchronization can be 
used to build a framework for network-wide synchronization 
across a sensor network. 



Appendix A 
Proof of Theorem|5] 

The ratio of the likelihood functions can be expressed as 
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The expectation of the ratio of the likelihood functions can 
now be calculated as 
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where it follows from ( fTTj ) that 

(M v(z) (h))- 2N = e-WHrt) 
Since the samples Zj are i.i.d, 
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With £(.) defined in the theorem, the proof is complete. 
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Proof of Theorem|7] 
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Continuing with the calculations 

T h (h k )= S (h k ) r /(p J fc + hfc)2 rfp, 
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It can be easily verified that ([84|> can be written as 



k /"+0O 

5(h fc ) = n / 



f(Zj\pjJ-hj] 
f{Zj\Pj) 



f(Z j \p j )dZ j 



Moreover, it can be noted that 

2 



{ f{Z f(Z-\l) 3) ) =^P[-2(^(ft+^)-^(P3))]x 
exp(2hjt] p (Zj)) 

and therefore 



+00 



fiHPi) 



f{Z 3 \ P j)dZ 3 =M-f{h 3 )K 



E [exp (2hjT)p(Zj))] . 

Then, since 

E [exp {2h jrip {Z 3 ))] = exp (4> p ( Pj + 2hj) - <p p (pj)) ( 85 > 
it can be easily seen that 



+00 



fjZjlpj + hj) 



f(Z J \p J )dZ J =M- 2 (h J )M Vp (2h J ) 



thus getting 



S(h k ) = ]jM- 2 (h j )M rip (2h j ) 
3=1 



Appendix C 
MSE Expressions for ML Estimators 



A. Gaussian Distribution 



If the likelihood for £ is Gaussian distributed ( fT~9] >, the MLE 
is given by pTj ). Since the variance of the readings Uj is cr| 
and the MLE ( pl~| is unbiased, it is straightforward to see that 



MSE Uml = 



N ' 



and similarly for ?/>ml- Given ( fT~8^ > and Proposition [l] it can be 
concluded that 



MSE ( 0v 



47V 



(86) 



B. Exponential Distribution 

If the likelihood for £ is exponential distributed pTj ), the 
MLE is given by ( [33) ). Through simple algebra it can be seen 
that C/(i) is exponentially distributed with parameter A ? = 

X^N, so that Var ^<|ml) = X 2' N 2 ■ It can be noticed that |ml 
is a biased estimator for £, with bias fo^ML = xTv- Similarly, 



Var (v>ml) = x^ap and Hml = j^n- Therefore, given ^ 
and Proposition|f] it can be concluded that 



MSEKl)=^(^ 




(87) 
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